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ABSTRACT 

O A modified non-linear time series analysis technique, which computes the cor- 
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relation dimension D 2 , is used to analyze the X-ray light curves of the black hole 



5_i ■ system GRS 1915+105 in all twelve temporal classes. For four of these temporal 

classes D 2 saturates to ~ 4 — 5 which indicates that the underlying dynamical 
mechanism is a low dimensional chaotic system. Of the other eight classes, three 
show stochastic behavior while five show deviation from randomness. The light 
curves for four classes which depict chaotic behavior have the smallest ratio of 
■ the expected Poisson noise to the variability (< 0.05) while those for the three 

classes which depict stochastic behavior is the highest (> 0.2). This suggests that 
the temporal behavior of the black hole system is governed by a low dimensional 
chaotic system, whose nature is detectable only when the Poisson fluctuations 



are much smaller than the variability. 
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1. Introduction 

Black hole X-ray binaries are variable on a wide range of timescales ranging from months 
to milli-seconds. A detailed analysis of their temporal variability is crucial to the understand- 
ing of the geometry and structure of these high energy sources. Such studies may eventually 
be used to test the relativistic nature of these sources and to understand the physics of 
the accretion process. The variability in different energy bands is generally quantified by 
computing the power spectrum which is the amplitude squared of the Fourier transform. 
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The power spectra give information about the characteristic frequencies of the system which 
show up as either breaks or as near Gaussian peaks, i.e Quasi-Periodic Oscillation (QPO) 
in the spectra (e.g. Belloni et al. (2001); Tomsick & Kaaret (2001); Rodriguez et al. 
(2002)). The shape of the power spectra, combined with the observed frequency dependent 
time lags between different energy bands, have put constraints on the radiative mechanisms 
and geometry of emitting regions (e.g. Nowak et al. (1999); Misra (2000); Cui (1999); 
Poutanen & Fabian (1999); Chakrabarti & Manickam (2000); Nobili et al. (2001)) 

These results are based on the response of the system to temporal variations whose 
origin is not clear. Important insight into the origin can be obtained by the detection 
and quantification of the possible non-linear behavior of the fluctuations. For example, 
the presence of stochastic fluctuations would favor X-ray variations driven by variations of 
some external parameters (like the mass accretion rate), or the possibility that active flares 
occur randomly. On the other hand, if the fluctuations can be described as a deterministic 
chaotic system, then inner disk instability or coherent flaring activity models will be the 
likely origin. A quantitative description of the temporal behavior can also be compared with 
time dependent numerical simulations of the accretion process and will help examine the 
physical relevance of these simulations. 

The non-Gaussian and non-zero skewness values of the temporal variation of the black 
hole system Cygnus X-l suggested that the variations are non-linear in nature (Thiel et 
al. 2001; Timmer et al. 2000; Maccarone & Coppi 2002). More rigorous tests were 
applied to the AGN ArK 564 (Gliozzi et al. 2002) which also suggested non-linear behavior. 
Nonlinear time series analysis(NLTS) seems to be the most convenient tool to check if the 
origin of the variability is chaotic, stochastic or a mixture of the two and has been adopted in 
several disciplines to study complex systems (e.g. human brain, weather) and predict their 
immediate future (Schreiber 1999). This technique has also been used earlier to analyze 
X-ray data of astrophysical sources. Based on a NLTS analysis of EXOSAT data, Voges et 
al. (1987) claimed that the X-ray Pulsar Her X-l was a low dimensional chaotic system. 
However, Norris & Matilsky (1989) pointed out problems with that analysis since the source 
has a strong periodicity and the data analyzed had low signal to noise ratio. Lehto et al. 
(1993) used the NLTS technique to analyze EXOSAT light curves of several AGN, and 
found that only one, NGC 4051, showed signs of low dimensional chaos. A similar analysis 
on the noise filtered Tenma satellite data of Cyg X-l, suggested that the source may be a low 
dimensional chaotic system with large intrinsic noise (Unno et al. 1990). These analysis 
were hampered by small number of data points (^ 1000) in the light curve and/or noise. 
Hence, the reported detection of low dimension chaos was only possible by rather subjective 
comparison of the results of the data analysis with those from simulated data of chaotic 
systems with noise. 
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The Galactic micro quasar GRS 1915+105 is a highly variable black hole system. It 
shows a wide range of variability (Chen et al. 1997; Paul et al. 1997; Belloni et al. 1997a) 
which required Belloni et al. (2000) to classify the behavior in no less than twelve temporal 
classes. In this work, our motivation is to determine the temporal property of this source by 
using a modified nonlinear time series analysis for each of these twelve classes. The different 
kinds of variability and its brightness ( the average RXTE PCA count rate ranges from 
5000 — 32000 counts/s) makes this source an ideal one to detect chaotic behavior. 

In the next section we describe the technique used to determine the Correlation dimen- 
sion. The results of the analysis are presented in §3, while in §4 the work is summarized and 
discussed. 



2. The Non-Linear time series analysis 

The algorithm normally employed in this analysis (Grassberger & Procaccia 1983) aims 
at creating an artificial or pseudo space of dimension M with delay vectors constructed by 
splitting a scalar time series s(t) with delay time r as 

x(t) = [s(t), s(t + r), , s(t + (M - l)r)] (1) 

The correlation sum or the correlation function is the average number of data points 
within a distance R from a data point, 

1 N N 

Cm{r) " J™ wn^T) £ £ m - - f ^ (2) 

where Xj is the position vector of a point belonging to the attractor in the M-dimensional 
space, N is the number of reconstructed vectors and H is the Heaviside step function. The 
fractional dimension D 2 (M) is defined as 

D 2 = lim dlogC M {R)/dlog{R) (3) 

— ^0 

and is essentially the scaling index of C M (R) variation with R. D 2 (M) can be used to 
differentiate between different temporal behavior since for an uncorrelated stochastic system, 
D 2 ~ M while for a chaotic system, D 2 (M) m constant for M greater than a certain 
dimension M max . 

For a finite duration light curve, there are two complications that hinder the successful 
computation of D 2 (M). First, for small values of R, C M (R) is of order unity and the 
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result there would be dominated by Poisson noise. Second, for large values of R, Cm(R) 
will saturate to the total number of data points. Usually, these two effects are avoided in 
the logC^-R) versus log/? plot and the slope D 2 is obtained from the linear part of the 
curve. However, such an exercise is subjective especially for high dimensions. Here, we use 
a numerical scheme to compute D 2 , which takes into account the above effects and at the 
same time optimizes the maximum use of the available data. The details of the method 
and several tests of its validity will be presented elsewhere ( Misra et al. in preparation). 
Briefly, the technique involves converting the original light curve to a uniform deviate, and 
to redefine the correlation function Cm(R) as the average number of data points within a 
M-cube (instead of a M-sphere) of length R around a data point. Only those M-cubes are 
considered which are within the embedding space, ensuring that there are no edge effects due 
to limited data points. This imposes a maximum value of R < R max for which Cm{R) can 
be computed. To avoid the Poisson noise dominated region, only results from R greater than 
a Rmin are taken into consideration such that the average C(R m i n ) > 1 where the Poisson 
noise would approximately be l/\fW c . Typically Cm(R) is computed for ten different values 
of R between R m i n and R ma x and the logarithmic slope for each point is computed and the 
average is taken to be D 2 (M). The error on D 2 (M) is estimated to be the mean standard 
deviation around this average. It should be noted that there often exists a critical M cr for 
which R max ~ Rmin and no significant result can then be obtained for M > M cr . Figure 1 
(a) shows the D 2 (M) curve for a time series generated from random numbers and for the 
well known analytical low dimensional chaotic system, the Lorenz system. The total number 
of data points used to generate both curves is 30000 and the number of random centers used 
is N c = 2000. As expected the D 2 plot for the random data is consistent with the D 2 = M 
curve, while the plot for the Lorenz system shows significant deviation and saturates at 
M ps 3 to a D 2 ps 2, which is close to the known value of 2.04. The random data and the 
low dimensional chaotic system can clearly be distinguished in this scheme. 



3. Results 

The temporal property of GRS 1915+105 have been classified into twelve different classes 
by Belloni et al. (2000), who also present the observational dates and identification num- 
ber of the RXTE data they had used to make the classification. Here, we have chosen a 
representative data for each class and extracted a few continuous data streams (ps 3000 sec 
long) from it. The observational IDs of the data used in this work are tabulated in Table 1. 
The light curves were generated with a resolution of 0.1 seconds resulting in ps 30000 data 
points for each of them and ps 1500 counts per bin. Light curves with finer time resolution 
are Poisson noise dominated, while larger binning gives too few data points. 
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In general, D 2 (M) is proportional to r when r is small and saturates (i.e. it is nearly- 
invariant) for r greater than a critical value and it is this saturated value which is the correct 
estimate of D 2 (M). As an example, the D 2 (M) curves for different values of r are plotted in 
Figure 1 (b), where it can be seen that the curve is similar within error bars for r = 15, 25 
and 100 sec. For all the data analyzed here the critical r < 5 — 20 sec, and hence the saturated 
curve (typically for r ps 50 sec) is considered. It has been verified that the D 2 (M) curves for 
two separate light curves for the same class, are similar to within the error bars. This shows 
that as expected the temporal behavior of the system is more or less stationary for the same 
class. Hence such curves can be averaged to obtain a statistically more significant result. 

Figure 2 shows the D 2 (M) curves for seven temporal classes. For four classes (A, k, f3 and 
n) the curves show clear deviation from random behavior. For A and k there is saturation of 
D 2 5 for M > 8. For (3 and /i, the increase in D 2 is less than one when M increases from 8 
to 15. Thus these classes can be classified unambiguously as chaotic systems with correlation 
dimension less than 5 while the the behavior of the class <fi is identical to a stochastic light 
curve. The classes a and p show some deviation from stochastic behavior and hence this 
behavior, which is also seen in the classes 9, v and 8, is named "non stochastic" in this work. 
As discussed below, these classes may be inferred to be low dimensional chaotic systems 
based on comparison with results from simulated data of chaotic systems with additional 
noise. Similar comparisons were made to infer the chaotic behavior of Cyg X-l (Unno et 
al. 1990) and NGC 4051 (Lehto et al. 1993). We show in the last column of Table 1 
the classification of all the twelve classes into one of these three categories, namely chaotic, 
non-stochastic and stochastic. We have listed in Table 1 the average counts < S >, the 
root mean square variation, the expected Poisson noise < PN >= V < S >, and the ratio 
of the expected Poisson noise to the actual RMS value. It can be seen that there is a 
strong correlation between the inferred behavior of the system and the ratio of the expected 
Poisson noise to the rms values. This indicates that Poisson noise is affecting the analysis. To 
estimate the effect of Poisson noise, we consider the Lorenz system points Sz,(t) and rescale 
it by Slr = AS lit) + B. A light curve is then simulated using Slr from the corresponding 
Poisson noise distributions. The constants A and B were chosen such that the simulated 
light curve had the same average count and rms variation as the two extreme cases for the 
GRS1915+105 data for f3 and 7 classes. The results of the non linear time series analysis are 
shown in Figure 3, where it can be seen that even for the f3 like case where the ratio of the 
expected Poisson noise to rms variation is only 4%, the D 2 versus M curve saturates at a 
higher value than that of the original no noise data points. This implies that the correlation 
dimension of ~ 4 inferred from the analysis of the classes showing chaotic behavior (Figure 
2) is an overestimation due to the inherent Poisson noise in the data. For larger Poisson noise 
fractions, the curve no longer saturates and becomes qualitatively similar to that obtained 
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for the non-stochastic case. 



4. Discussion 

The saturation of the correlation dimension D 2 ~ 4 — 5 for four of the temporal classes 
clearly indicates that the underlying dynamic mechanism that governs the variability of the 
black hole system is a low dimensional chaotic one. As indicated by simulations of the 
Lorenz system with noise, the effect of Poisson noise in the data is to increase the D 2 values. 
Hence the real dimension of the system is probably smaller than D 2 ~ 4 — 5 that is obtained 
here. In fact it is possible that the the temporal behavior of the black hole system is always 
governed by a low dimensional chaotic system, but is undetectable when Poisson noise affects 
the analysis. 

Alternatively, there may be a stochastic component to the variability which dominates 
for certain temporal classes. The two scenarios may be distinguished and better quantitative 
estimates of the correlation dimension may be obtained by either appropriate noise filtering 
of the data and/or appropriate averaging of the different light curves. Much longer (?» 30000 
sec long) continuous data streams sampled at 1 second resolution, would decrease Poisson 
noise and hence provide better quantitative measure of D 2 . However, such long data streams 
are presently not available and merging non-continuous light curves, will require sophisticated 
gap filling techniques which might give rise to spurious results. 

The variability of GRS 1915+105 can be interpreted as being transitions between three 
spectral states (Belloni et al. 2000), one of which (the so called soft state) is a long term 
canonical state observed in other black hole systems like Cygnus X-l which do not show such 
high amplitude variability. It is attractive to identify these spectral states as fixed points 
which for GRS 1915+105 become unstable giving rise to the observed chaotic behavior 
which may also account for the ring like movement of the system in color-color space (Vilhu 
& Nevalainen 1998). The above hypothesis may be verified by future characterization of the 
chaos in GRS 1915+105. Note that GRS15+105 spends most of it's time in the \ class whose 
variability is similar to that observed in other black hole systems like Cygnus X-l. However, 
as shown in this work, Poisson noise effects the analysis for the x class and the D 2 (M) values 
reflect stochastic behavior. This may be the reason why earlier different non-linear analysis 
of Cygnus X-l data, while showing non-linearity (Timmer et al. 2000; Thiel et al. 2001) 
did not conclusively reveal chaotic behavior. 

The identification of the temporal behavior of the black hole system as a chaotic one, has 
opened a new window toward the understanding of the origin and nature of their variability. 



-7- 



The present analysis can be extended to characterize the chaotic behavior. Using the mini- 
mum required phase space dimension, the data can be projected into different 2 dimensional 
planes, which will reveal the structure of the attractor and help to identify any possible 
centers of instability in the system. Further, dynamical invariants like the full Lyapunov 
spectrum , multi-fractal dimensions etc. can be also be computed. Recently, Winters et al. 
(2003) have studied and quantified the chaotic flow in magneto-hydrodynamic simulations 
of the mass accretion processes that is believed to be happening in black hole systems. The 
measured chaos parameters like the largest Lyapunov exponent for such simulations can be 
compared with that obtained from the light curve of black hole systems to validate such 
simulations and enhance our understanding of these systems. Note that such analysis can 
practically be applied only after the identification of the minimum phase space dimension 
which in turn usually requires the computation of D 2 (M). 

GA and KPH acknowledge the hospitality and the facilities in IUCAA. BM thanks the 
Academy of Finland grant 80750 for support. 
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Table 1. Observation Table 



Obs. I.D. 


Class 


< s > 


rms 


<PN> 


< PN > /rms 


Behavior 


10408-01-10-00 


P 


1917 


1016 


43.8 


0.04 


C 


20402-01-37-01 


A 


1493 


1015 


38.6 


0.04 


c 


20402-01-33-00 


K 


1311 


800 


36.2 


0.04 


c 


10408-01-08-00 


n 


3026 


999 


55 


0.06 


c 


20402-01-45-02 


e 


1740 


678 


41.7 


0.06 


NS 


10408-01-40-00 


V 


1360 


462 


36.9 


0.08 


NS 


20402-01-03-00 


p 


1258 


440 


35.5 


0.08 


NS 


20187-02-01-00 


a 


582 


244 


24.1 


0.10 


NS 


10408-01-17-00 


5 


1397 


377 


37.4 


0.10 


NS 


20402-01-56-00 


7 


1848 


185 


43.0 


0.23 


S 


10408-01-22-00 


X 


981 


118 


31.3 


0.27 


s 


10408-01-12-00 





1073 


118 


32.7 


0.28 


s 
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Fig. I. — (a) The D 2 versus M curve for random points (circles) and for Lorenz system 
(squares). For both the curves the number of points used are 30000 and the number of 
centers used in the computation is 2000. The straight line is the D 2 = M case which is the 
expected result for random variation, (b) The D 2 versus M curves for GRS 1915+ 105 data 
obtained during the class k for three different values of the delay time r = 15 s (triangles), 
= 25 s (squares) and = 100 s (circles). 




Fig. 2. — The D 2 versus M curves for GRS1915+105 data obtained during seven temporal 
classes. The classes re, /x, (3 and A exhibit chaotic behavior. The classes is stochastic, while 
a and p show some departure from stochastic behavior. 
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Fig. 3. — The effect of Poisson noise on the analysis is studied using simulated rescaled 
Lorenz system data with Poisson noise. The D 2 versus M curve for Lorenz system is shown 
by filled circles, while the squares (triangles) are curves for rescaled Lorenz system data 
corresponding to the expected Poisson noise in the (3 (7) class of GRS 1915+105 



